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A method is presented for the optimization of one-body and inhomogeneous two-body terms in 
correlated electronic wave functions of Jastrow-Slater type. The most general form of inhomogeneous 
correlation term which is compatible with crystal symmetry is used and the energy is minimized 
with respect to all parameters using a rapidly convergent iterative approach, based on Monte Carlo 
sampling of the energy and fitting of energy fluctuations. The energy minimization is performed 
exactly within statistical sampling error for the energy derivatives and the resulting one- and two- 
body terms of the wave function are found to be well-determined. The largest calculations performed 
require the optimization of over 3000 parameters. The inhomogeneous two-electron correlation terms 
are calculated for diamond and rhombohedral graphite. The optimal terms in diamond are found 
to be approximately homogeneous and isotropic over all ranges of electron separation, but exhibit 
some inhomogeneity at short- and intermediate-range, whereas those in graphite are found to be 
homogeneous at short-range, but inhomogeneous and anisotropic at intermediate- and long-range 
electron separation. 

PACS numbers: 71.15.-m, 71.15.Nc, 02.70.Uu 
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I. INTRODUCTION 

An accurate description of electron correlation is one 
of the central issues in modern electronic structure the- 
ory. In principle, this involves solving the many-electron 
Schrodinger equation, which for a system of N inter- 
acting electrons is an inseparable 37V-dimensional prob- 
lem. Methods based on mean-field approximations of 
the electron interaction, .such as Hartree-Fock (HF) or 
density-functional theory,El (DFT) reduce this to TV inde- 
pendent problems, which we may solve to approximate 
many physical properties with reasonable accuracy. How- 
ever, a more complete description of electron correlation 
is often required to study delicate phenomena such as 
long-range electron correlation and van der Waals inter- 
actions. 

Quantum Monte Carlo (QMC) methods provide an im- 
portant .approach for solving the full 3N dimensional 
problemE One such method, variational Monte Carlo 
(VMC), allows us to estimate expectation values for a 
given trial wave function. Ideally, this wave function is 
an eigenstate of the many-body Hamiltonian. In prac- 
tice, a parameterized form is used, which approximates 
the exact eigenstate. The accuracy of these calculations 
is entirely dependent on the trial wave function, however, 
and so the development of accurate wave functions is vi- 
tal both in the accurate estimation of physical properties 
and in understanding how certain physical phenomena 
may be simply represented in the wave function. 

A widely used trial correJated wave function is the 
Jastrow-Slaterta or Feenbergcl form, 



#(ri, ...,r N ) = exp 



i<3 



D.(l) 



Here I? is a Slater determinant of single-particle orbitals 
and interparticlc correlation is introduced with the two- 
body term u in the Jastrow factor. The one-body term 
X could in principle be absorbed into the single-particle 
orbitals of the determinant, but it may be convenient to 
retain it explicitly in the Jastrow factor. In practice, the 
orbitals in the Slater determinant are often determined 
from a HF or DFT calculation.0 

In this paper, we will focus on the optimization of 
Jastrow-Slater wave functions in the context of the elec- 
tronic structure of periodic solids. We will apply wave 
function optimization to examine some consequences 
which emerge from a complete treatment of inhomogene- 
ity in the two-body correlation term for diamond and 
rhombohedral graphite. 

The general form of wave function in Eq. [l] has been 
used as the starting point of several methods, includ- 
ing Fermi hypernetted chainO (FHNC) and VMC0 cal- 
culations. The traditional approach is to use a varia- 
tional principle on the energy (or, in the "variance min- 
imization" method^ on the fluctuations of the energy) 
to define a best approximation to the true eigenstate 
within the variational freedom allowed by the wave func- 
tion ansatz. Expectation values of the energy and other 
quantities are calculated from the trial wave function, 
approximately in the FHNC approach, exactly in VMC 
(within statistical error of the sampling). While the VMC 
method has been used in calculations of periodic solids, 
in practice most calculations have included only homo- 
geneous two-body terms in the Jastrow factom and the 
optimization of wavefunctions with very large numbers of 
parameters has remained problematical. To our knowl- 
edge, the FHNC approach has not been applied in fully 
three-dimensional electronic structure calculations. 
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Wave functions determined by the VMC approach are 
often used as guiding, trial functions for the diffusion 
Monte Carlo (DMC) method.Q In the DMC context, the 
accuracy of the wavefunction affects the numerical effi- 
ciency of energy calculations and the accuracy of other 
physical quantities. When DMCpis used in conjunction 
with non-local pseudopotentialsp as is commonly the 
case in chemical applications, an accurate trial wavefunc- 
tion is essential for an accurate calculation of the ground 
state energy. 

We present here, for the first time, a numerically ro- 
bust, rapidly convergent iterative method which mini- 
mizes the variational energy with respect to a very gen- 
eral inhomogeneous form of the Jastrow factor, includ- 
ing all one- and two-body terms compatible with crys- 
tal symmetry. The method remains numerically well- 
conditioned, even for large systems (of the order of sev- 
eral hundred electrons), where there are more than 3000 
independent variational parameters in the Jastrow fac- 
tor. Within acceptable computational demands on the 
Monte Carlo sampling, the optimal values of parameters 
in the wave function are found to be well determined, 
even when their contribution to the total energy is ex- 
tremely small. 

The method places no restrictions on the functional 
form of the anti-symmetric (determinantal) part of the 
many-body wave function and can be used in conjunc- 
tion with related methods recently developed for en- 
ergy minimization with respect to all orbitals in the 
determinantta or with respect te.rGonfiguration weights in 
a multi-determinant function.! 10 !! 11 ! Taken in conjunction 
with these methods, the approach presented here com- 
pletes the solution to the problem of energy minimiza- 
tion with respect to the most general variational terms 
in wavefunctions of the Jastrow-Slater type, as well as 
in wave functions where the Jastrow factor multiplies 
a multi-determinant function. Although many specific 
details of the work here refer to periodic systems, the 
basic method for Jastrow factor optimization could also 
be used in similar calculations of molecular, atomic, or 
nuclear structure. 

Other general methods exist to achieve energy mini- 
mization with respect to parameters in many-body wave 
functions. When only one or two parameters are op- 
timized, it is possible to perform a systematic search of 
parameter space, as McMillan did ip-his pioneering VMC 
study of the properties of liquid He.lI3 The stochastic gra- 
dient approximation (SGA) used by Harju et allEi applies 
control theory to determine iterative corrections to the 
wave function parameters. Lin et al.t3 explicitly com- 
puted analytical derivatives of the energy with respect to 
variational parameters and used these to optimize their 
wave function with a Newton-style method. However, 
none of these methods have as yet been applied in sys- 
tems where a very large number of parameters are to be 
optimized. 

In recent years, the variance minimization method of 
Umrigar et al.p which optimizes the wave function by re- 



ducing the magnitude of the variance of the local energy, 
has been used much more widely than energy minimiza- 
tion in determining optimal parameters for variational 
wave functions. Although not strictly equivalent to en- 
ergy minimization for a non-exact wave function, in prac- 
tice this method has been very successful in improving to- 
tal energies, providing extremely accurate wave functions 
in certain atomic systems.Q Unfortunately, variance min- 
imization is often subject to confinement to local minima 
and requires much human interaction and experience for 
successful implementation when large numbers of param- 
eters are to be optimized. 

Variance minimization may be thought of as fitting en- 
ergy fluctuations to a constant, and attempting to reduce 
the cost function involved in this fit by direct variation of 
the wave function parameters. The standard approach-is 
to use a least-squares fit of these energy fluctuations x3 

Our method of energy minimization involves the fitting 
of energy fluctuations to a given (non-constant) func- 
tional form, which is a linear combination of operators 
associated with variations of the wave function parame- 
ters (see Sec. III). This fitting allows us to determine a 
"predictor" , which links changes in the wave function to 
changes in the fitted energy fluctuations. This predictor 
iteratively guides our method to a self-consistent solu- 
tion, where the fitted energy fluctuations are zero. The 
derivatives of the true many-body energy with respect to 
all the parameters in the wave function are then also zero 
(within statistical sampling error) for this final solution. 

Our predictor is closely related to the random phase 
approximation (RPA), introduced by Bohm and PinesEj 
and recejuthi discussed in the context of inhomogeneous 
systems. Ej'EZI However, as long as the predictor is suffi- 
ciently accurate to ensure stable convergence of the iter- 
ations to the self-consistent solution, its exact form does 
not affect the final solution. Our method allows us to 
surpass the approximations of the RPA and produce ex- 
plicit trial wave functions of unprecedented accuracy for 
electrons in periodic solids. 

We will use the optimized wave functions to study the 
effects of charge inhomogeneity on the correlation factors 
in diamond, as a prototype of strongly bonded insulating 
systems, and rhombohedral graphite, as a prototype of 
highly anisotropic, inhomogeneous solids. Although ear- 
lier studiesErES of these and related systems would indi- 
cate that homogeneous two-body correlation factors gain 
a large fraction of the correlation energy in solids and 
that inhomogeneous correlation factors are unlikely to 
lower the total variational energy greatly, substantial dif- 
ferences in correlation factors often give rise to relatively 
little change in energy. This is particularly true when 
one is interested in long-range correlation, which is ener- 
getically very delicate. 

The rest of this paper is organized as follows: In Sec. 
U, we present the detailed numerical form of variational 
wave functions to be used in these calculations. The ap- 
proach of fitting energy fluctuations and its use in guiding 
the iterative solution of the energy minimization problem 
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is discussed in Sec. HI. In Sec. IV, we present some re- 



sults on the application of the method to the periodic 
solids, diamond and rhombohedral graphite, and exam- 
ine the effects of charge density inhomogeneity on the 
correlation factors of these systems. In Sec. [v] we dis- 
cuss the results and some computational details of the 
method, illustrating the justification for certain aspects 



of our approach with some tests. Finally, in Sec. VI 
present the overall conclusions of this study. 



we 



II. FORM OF THE WAVE FUNCTION 

The Jastrow two-body correlation factor u(r, r') in Eq. 
[I] can in principle always be expressed as: 



a/3 



(2) 



where f a form a complete set of functions and u a p are 
expansion coefficients. (We use * to indicate complex 
conjugation throughout this paper). Similarly, the one- 
body function x( r ) may also be expanded in a basis set 
of complete functions g 7 , as: 



(3) 



Apart from-an additional term to handle the electron- 
electron cusp,tj we will express the Jastrow factor in the 
general form given by Eqs. || and [| Values of the parame- 
ters {u a /3, Xj} will then be determined so that the energy 
of the system is stationary with respect to all variations 
{du Q/3 ,dx 7 }- 



A. The electron-electron cusp 

Due to the divergence of the Coulomb interaction be- 
tween two electrons, the correct two-body correlation 
term u(r, r') has a cusp where r — > r', leading to slow 
convergence of-any expansion in smooth functions of the 
form in Eq. |^.Ey For this reason, it is numerically conve- 
nient to re-write the two-body function u in the form: 



M(r,r')+M sr (|r-r'|), 



where u ST (r) is a short-ranged (homogeneous) function 
which has the correct electron-electron cusp as r — > 
and it(r,r') is now a smooth, cuspless function. We use 
a form of the short-ranged function u ST which is gener- 
ated from a numerical solution of the electron-electron 
scattering problem. A discussion of the generation of u ST 
is provided in Appendices A and C of Ref. We de- 
fine u sr = — In J sr , where J sr is defined within Appendix 
C of Ref. Ell The expansion of the remaining function 
u(r, r') in the form given in Eq. ^| then converges much 
more rapidly than that of the original function, for any 
set of smooth functions f a . 



In this paper, we generate the short-range function 
u sr in a spin-dependent form, to maintain the cusp 
conditions .113 The cuspless function u(r,r') used in our 
work is independent of electron spin, although we expect 
that a spin-dependent form is as easily optimized. 



B. Separation of one- and two-body terms 

Summing u(r, r') over all electron pairs in the basis f a 
leads to 



(5) 



i<j 



a.0 



This may be thought of as a two-body expansion of the 
correlation term. However, it is important to realize that 
any separation of "one-body" and "two-body" terms in 
the Jastrow factor is somewhat arbitrary. Removal of 
terms where i = j, as in Eq. |^, is not sufficient to de- 
couple one- and two-body terms completely. To see this, 
consider the transformation of each of the basis functions 
obtained by subtracting a constant, f a (r) — f a (r) — c a . 
The set remains complete and the function 

u(r, r') = £[/ Q (r)* - c* a ]u a ^(r') - c p ] , (6) 

a/3 

may be interpreted as a "two-body" function in the basis 
set f' a . However, expanding this correlation factor over 
all electron pairs, we find that, in terms of the original 
basis f ai an additional one-body contribution appears: 



a0 



N -1 



^^[u a[j f a {v i )*ci3 + u a f i c a ff S (v t )]{l) 

af3 i 



constant 



We may regard this as a transformation of the one-body 
function in Eq. [|: x( r ) ~ > x( r ) + X°( r )> where the addi- 
tional term comes from the second line of Eq. ^ 



(4) x°(r) 



N — 1 



(8) 



E 











fair) 



To make a definite numerical separation of our one- 
and two-body expansions, we need to consider an appro- 
priate choice of the arbitrary constants c a , noting that 
u(r, r') is intended primarily to affect correlation prop- 
erties (i.e., two-body properties) of the system, leaving 
single-particle properties unchanged. For example, the 
mean- field methods that produce D [HF, DFT under the 
local density approximationE3 (LDA), etc.] normally give 
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very accurate single-particle densities, which can be al- 
tered substantially by inclusion of an arbitrary function 
u(r, r'), in 'f (Eq. |l|). Any substantial change in the den- 
sity from the HF solution is likely to be energetically very 
costly, so that ideally we would like to decouple changes 
in it(r, r') from changes in the density. Necessary changes 
in the density may be allowed for by optimization of the 
explicit one-body term \ in the Jastrow factor (Eq. ||) or 
by the methods of Ref. ||. 

Therefore, we would like to choose the constants c a 
such that the average of any one-body operator (such 
as the density) for W remains stationary with respect to 
variations in the coefficients u a p, at least in the absence 
of interparticle correlation. If we define the one-body 
operator 

i 

for the many-body configuration R = (ri ...rjv), then 
its expectation value is 

(p 7 ) = (vE-KI*) = J / 7 (r)*p(r) dr , 

where p(r) is the single particle density. Since the / 7 
are fixed functions, variations in the expectation values 
(p 7 ) correspond to variations in p(r). The derivative of 
(p 7 ) with respect to variations of the u a /3 in Eq. ^| is (see 
Appendix [X|) 

\jt^ = (M R ) - (pj)\ £[A( r *r - cWfa) - <*]> 

k 

In the absence of interparticle correlation, 

r (f;f*)(M nk=i 
(f^ruv^Mr,)) = { (/*//?></<*>* if k = j 

{ (fy)*(Ur(U) otherwise 

The second summation excludes i — j, so averages of 
triple products never arise. Thus, in the absence of cor- 
relation, the derivative becomes 

I^W ~ (9) 

2 dUa/3 

£i*y (IM^y ~ (/ 7 >1[/a(r<)* " c^iMr,) - cp) 

We can guarantee that the right hand side of Eq. ^ is 
zero if c a = (f a ), for all a. In other words, one-body 
expectation values remain approximately unaffected by 
the presence of the correlation factor w(r,r'), provided 
we expand u(r, r') in a basis of "fluctuation functions", 
fa(r) = f a {r) — (f a ). Equivalently, we may retain the 
original basis f a and form a one-body term x° from Eq. 



H with c a = (f a )- This may then be inserted into the 
wave function of Eq. [j] and varied as the parameters u a p 
are varied. 

Correlation effects are of course present in the actual 
wave function used. However, we find that Eq..Kj remains 
approximately true, as previously observed.E3E3 For the 
energy minimization problem, we find that mixed deriva- 
tives of the energy d 2 E / du a pdx-i are approximately zero 
when c a — (/ Q ). This gives the numerical advantage 
that minimization of the energy with respect to u a p ap- 
proximately decouples from minimization with respect 
to X 7 . We note that this approximate decoupling holds, 
even if the relation c a — (f a ) is not exactly true. Thus 
we may use c Q = (D\f a \D), (i.e., LDA or HF aver- 
ages of f a ), in place of (\E , |/ Q |^I'), while still maintain- 
ing the numerical advantages of approximately satisfying 
d 2 E/du a pd\ 1 = 0. 

C. Fourier Expansion 

In the context of periodic systems, it is natural to ex- 
pand the correlation function u(r, r') as a Fourier series, 
where the basis functions are / q = exp[iq- r], for each 
wave vector q. We note that 

p q (R) = £ expHq • rj = £ / q (r,)* (10) 

i i 

is the Fourier coefficient of the instantaneous charge den- 
sity, given the electron configuration R. (In atomic units 
the electronic charge e = 1). Also, summing over pairs 
leads to a quadratic product of Fourier coefficients, with 
some modification to remove terms with i = j: 

EAWA'fo) = ^(PqPtfW (11) 

i<j 

In order to approximately remove the effect of the two- 
body terms on the single-particle density, following Sec. 
15 B| , we subtract appropriate constants from each basis 
function to produce a new basis, the collective "charge 
fluctuation" coordinates, 

E U r *y = = pi (pi) = E fi^ </i>* ( 12 ) 

i i 

These provide a suitable expansion of the two-body corre- 
lation factor that approximately preserves single-particle 
densities. 

A correlation factor using such coordinates was first 
suggested by Bohm and PinesEj for the homogeneous 
electron gas, and has recently been discussed by Malat- 
esta et oIIlB and Gaudoin et alHi in the context of inho- 
mogeneous systems. In homogeneous systems the expec- 
tation value (p q ) of each charge density Fourier coefficient 
is zero, and so charge density fluctuations are simply p q . 
For inhomogeneous systems, in general (p q ) ^ when 
q = G, a reciprocal lattice vector. 
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The alternative to using fluctuation coordinates Ap q 
is to incorporate the equivalent one-body term in the 
Jastrow factor, which in periodic systems is of the form 



,x° G Ph 



G 



with the coefficients coming from the two-body term 



v° 

AG 



N 



N 



— ^2 u gg>(pg> 



(13) 



G' 



as discussed by Malatesta et al. El and Gaudoin et alH 

The properties of the Fourier basis / q and the corre- 
lation factor u lead to some convenient symmetry prop- 



erties for the coefficients 



The complex conjugate 



M qq' = u -q-q'' J ust as /q = /-q' Tne exchange sym- 
metry of u, i.e., u(r,r') = w(r',r) implies that u qq / = 
w_ q / _ q . And, if u possesses inversion symmetry, i.e. 
w(r,r') = u{— r, — r'), then each u qq ' is a real number. 
In periodic systems, 

u(r + L,r' + L) = u(r,r') 

for any Bravais lattice vector L. This implies that all 
Fourier coefficients u qq ' are zero unless q q' = G, a 
reciprocal lattice vector. Thus, translation symmetry 
greatly reduces the number of variational parameters in 
the two-body terms. 

We arrange the wave function in the form 

f = J sl J ih J lh D, 

isolating the short range component of the Jastrow factor 
as 



J sr = exp 



E V 

|G|<G C 



G 



i<j 



r {Tij ) 



The one-body term here is derived from the short range 
correlation factor u sr of Eq. |], as Xg r (G) = u sr (G)(pG), 
where M sr (G) is the Fourier transform of u sr for the re- 
ciprocal lattice vector G. The prefactor of (N — 1)/N, 
which should be present from Eq. [l3|, approaches unity 
for large systems and so is neglected. For computational 
efficiency we leave this one-body term in its Fourier space 
representation and G c is a cut-off chosen for the Fourier 
sum such that it is converged within a required accuracy. 

The remaining inhomogeneous part of the two-body 
Jastrow factor is expanded using charge density fluctua- 
tion coordinates: 



J ih = exp 



} , } , "q+G q+G'-Pq+C q G 
q GG' 



, (14) 



where P q +c q +c = (A/3 q+G Apa±G')[i#j]) usin S tnc no- 
tation as defined in Eq. [Ll| and the definition of 
A/9 q _|_G m EQ- In practice, this double sum is trun- 
cated by using vectors q + G of magnitude less than a 
suitably chosen cut-off k c . 



We also allow for one-body optimization through the 
use of the explicit one-body Jastrow factor. This one- 
body Jastrow factor is also expanded in fluctuation co- 
ordinates, 



Jib = exp 



E 

|G|<G C 



(15) 



since including the constant average values (pg) merely 
adjusts the normalization of the wave function. 

The coefficients it q +G q +c and \g, defined in Eqs. 
|l4| and are the final variational parameters of our 
wave function. The remaining sections of this paper de- 
scribe the method we use to optimize these parameters 
such that the total energy of a particular electronic sys- 
tem is stationary. A typical calculation presented below 
involves the simultaneous optimization of over 3000 pa- 
rameters. 



III. ENERGY MINIMIZATION 

We wish to optimize the wave function 
* = *(a) , 

where a = {a m } a vector of parameters, by solving the 
Euler-Lagrange equations, 



d(H) 
da m 



= for all 



(16) 



We note that the Hamiltonian for the system is 

n = \ E v " + E W r + E v tra) > 



i<j 



where each sum is over the electrons in the system. 

Solving Eq. [l6| is equivalent (see Appendix |a]) to solv- 
ing the following system of equations 



(m AC m ) = for all m , 



(18) 



where, given a many-body configuration R = {r^}, we 
define: AA(R) = A(R) — (A), for any operator A; and 
the local values of the operators O m as, 



O m (R) 







ln*(R) (19) 

-P qq '(R) for a rn = M qq / , 
Ap G (R)* for a rn = \g ■ 



We shall refer to the local value of the Hamiltonian op- 
erator as the "local energy", 



E(R) = ™W 



(20) 



We approach the problem of solving the Euler- 
Lagrange equations (Eq. |l6| ) indirectly, by considering 
systematic fluctuations of the energy for a given trial 
wave function ^>(a). 
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A. Systematic Energy Fluctuations 

Consider fitting the local energy E(R) to the func- 
tional form 



E 



y m o m (R) 



(21) 



in the least-squares sense, where {O m } is the set of func- 
tions with which we fit the energy, and { Vm} is the vector 
of fitting coefficients. The least-squares problem reduces 
to minimizing the integral 



which is equivalent (see Appendix |b|) to solving the linear 
system 

53 V m (AO m AO n ) = (AEAO n ) for all n. (22) 



We recognize immediately that if the functions O m are 
those functions associated with variations of the wave 
function parameters a m (Eq. |l9| ), then the right-hand 
side of Eq. 22 is the vector of Eulcr-Lagrange derivatives 
in Eq. [l8| Therefore, the Euler-Lagrange equations (Eq. 
|l6| ) are solved if all the fitted coefficients V m are zero.c3 

As an illustrative example, consider & n eigenstate 
of TL. Indeed, the local energy is a constant, indepen- 
dent of R, and so we would find that each of the fitted 
coefficients V m is zero. Therefore, the Euler-Lagrange 
derivatives (AHAG m ) are all zero, and so the energy 
must be stationary with respect to variations in ^Sq, as 
we would expect for an eigenstate. 

For the trial wave function ty(a), no choice of a gives 
an exact eigenstate of the Hamiltonian. However, for a 
particular choice of parameters, the absence of system- 
atic variations of the energy (i.e., variations correlated 
with the variations of the functions O m ) ensures that 
the fitting coefficients V m are all zero and that the aver- 
age energy is stationary with respect to all variations of 
the parameters a.. Within the parametric freedom of the 
trial wave function '5 , this is our best approximation to 
an eigenstate. 



B. Iterative Procedure 

We now describe a procedure which aims, by appro- 
priate choice of the parameters a, to set the fitted co- 
efficients V m of the total energy E to zero. As defined 
in Eq. ^2], these V m depend on the wave function 'I' (ex) 
and so are functions of the parameters a. However, the 
functional dependence of V m on a is not available in its 
exact analytic form and we are unable to solve directly 
the system V m (ot) = 0, i.e., to find the root a which 
will guarantee the solution to the corresponding Euler- 
Lagrange equations. 



Instead, using the wave function ^(a°), for a particu- 
lar choice of the parameters a , we construct a predictor 
function V^ot.] oft) , which approximates the unknown 
function V m (a) for general values of a. More precisely, 
we construct V^(a;a°) so that 



V^( a °;a°) = V m ( ( x ) 
(i.e., V^j is exact when a = or) and, 



(23) 



(24) 



for all relevant values of a. 

To determine this predictor, we use the specific form 
of the Hamiltonian (Eq. [It]) and trial wave function ^ 
(Sec. ||), and we partition the local energy E(R) into a 
sum of contributions, 



25(R) = £eW(R) , 



(25) 



where the come from various terms in the kinetic and 
potential energy (see below). Each contribution eW(R) 
is approximated with the functional form 



v%O m (R) 



(26) 



For some terms, using the specific form of the local en- 
ergy for \I/(a), we can expand eW analytically as in Eq. 
p6| , enabling us to determine, exactly or approximately, 
the function Vm{oc). In this analytic form, Vm (a) is in- 
dependent of the choice of a and so remains equally 
valid for all a. 

Where analytic expressions are too complex to derive, 
we may approximate Vm (ct) by fitting to Eq. |2^. 
The fitting coefficients for e^ l \ are found by solving the 
analogue of Eq. p2|, 



where Vm — Vm (a ) is determined by using ty(a°) to 
evaluate the required expectation values. This produces 
the value of each coefficient at a . We may also fit the 
derivatives of with respect to a, to determine the 
linear dependence of each Vm(a). We may then approx- 
imate the function 



(AO m AO„) = (Ae^AOn) for all n, (27) 



.(0/ 



.(<) 



(a) 



t,«(a;a )=*,«(a ) + ^ 



dv 



(i) 



I 



da> 



(28) 
(a; - a?) , 



where the values of the derivatives are found by fitting 
de^/dai to Eq. H 

using ^>(a°). In evaluating the term 
dvm /dai, we consider only the explicit variation of the 
term Vm with on. We do not include the implicit varia- 
tion due to the dependence of the probability distribution 



7 



^(a)! 2 on a. In practice, the only term for which we 
need to fit de^ /dai is explicitly linear in the parameters 
a and so the linear expansion in Eq. ^8] is valid over a 
wide range of values of a. 

Just as the energy contributions (R) partition the 
local energy, we may regard the analytic and fitted coeffi- 

(i) 

cients Vm as an approximate partition of the local energy 
coefficients V m . We define this partition as follows: 

W m {a; a ) = S m (a) + T m {a; a°) , 

where the sum of analytically derived coefficients is 

S m (a)= £ 

analytic 

and the sum of numerically determined coefficients, eval- 
uated using ^(a ) and Eq. |27], is 



fitted 



We construct the predictor V^(a; a ) such that it satis- 
fies Eq. |2^, i.e., 

= V m (<x°)+W m (a;a°) - W m (a ;a°) (29) 



where Vmic* ) are found by solving Eq. 22 



We define our iterative approach to determining the 
parameters a, for which the coefficients V m (a) are zero, 
as follows: 



1. Given the set of parameters ot n , construct the wave 
function 4 r (cr n ). 

2. Evaluate the required expectation values in Eqs. |22| 
and |2^| using 4 r (ct n ), to find the fitting coefficients 
V m (a n ) and numerical functions Vm(a;a n ). 

3. If the total energy fitting coefficients are zero, i.e., 
V m (a n ) — for all m, then we are done, otherwise 
continue. 

4. Construct the predictor function V^Jat; a n ) in Eq. 
p9| using the fitted terms from Step. |2|, and find the 
solution ot n+1 to the system 



n+1. 



for all m, 



(30) 



using the Newton-Raphson methodo (see Ap- 
pendix^). Use this set of parameters a n+1 in Step 



Iterations continue until the total energy coefficients 
V m (a n ) tend to zero, and the values of the parame- 
ters a™ converge. Note that even though the predictor 
V^ n (a; a n ) is only an approximation to the exact function 
V m (a), Eq. 23 guarantees that at the converged solution 



In other words, the parameter set a solves the Euler- 
Lagrange equations for (H) exactly. Clearly, the larger 
the neighbourhood within which the approximate rela- 
tion in Eq. ||J holds, the faster this iterative procedure 
will converge. In the trivial case, if the exact analytic 
form of V m (ot) were known a priori, then we could just 
solve the Euler-Lagrange equations in one step using a 
suitable root-finding method. 

Rather than starting the procedure from an initial 
guess a = 0, we may begin at Step || using only the 
analytic terms in the predictor, since these are indepen- 
dent of the wave function and so do not require fitting. 
That is, we solve the system 

S m (at ) = for all m . 

The solution set a 1 is then used in Step [j]. 



C. Partitioning the local energy 

We note that the potential energy operators T4 x t and 
V in the Hamiltonian TL, defined in Eq. [l?], are multi- 
plicative, and therefore their contributions to E(R) are 
constant with respect to variations of the wave function 
parameters a. Variations of a affect only the contribu- 
tions of the differential kinetic energy operator. Thus, 
solving the Euler-Lagrange equations for (H) (Eq. |l6| ) 
amounts to adjusting the systematic fluctuations of the 
kinetic energy to cancel those in the potential energy ex- 
actly. 

If we extract a variational part <j>{a) of "J, such that 
$ = (fi^'i where <j> is dependent on the set of parameters 
a, and is independent of them. Then we may partition 
E(R) as follows: 

E(R) = e (1) (R) + e (2) (R) + e (3 >(R) 

where we define 



2 h *m 



.(2) 



ViHR) v^'(R) 
^ 0(R) ' *'(R) 



iV 



i—1 i—1 i<j 

Clearly, is constant with respect to variations of a. 
Further analysis of each of and is necessary to 
determine how each depends on a. However, for (f> ex- 
pressible in the form of the Jastrow factors in Sec. |J, 
i.e., 



N 



N 



= V^{a;oc)=V m (a) 



exp 
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where Oi are defined in Eq. E^, we see that: (i) 
is at most linear in a, but involves terms coming from 
so that it may be impossible to determine analytic 
expressions for the coefficients (a); (ii) e^ 1 ^ is at most 
quadratic in a, and involves only </>, for which we have 
an analytic expression, and therefore may derive analytic 

„(!)/ 



approximations to the coefficients v m ' (a) (see Sec. Ill D 



D. Analytic terms in the predictor 

The initial predictor S m (cx) used in the iterative 
method involves only analytic terms. We determine these 
by direct expansion of particular local energy contribu- 
tions, given the analytic form of the wave function. We 
now consider the contributions of the one- and two-body 
terms individually. 



One-body Jastrow Factor 



Replacing 4>(a) in Sec. [II C with Jib(x), we derive an 
analytic expression for the energy contribution e' 1 ' (see 
Appendix |D|). From this we may extract those coeffi- 
cients Wq" 1 of the functions Cg = A/?g : 



By assumption, the one-body contribution of e^ 3 -* is zero. 
That is, the mean field methods used to calculate the 
determinant D should remove (approximately) all sys- 
tematic one-body fluctuations in the local energy, and 
the use of "fluctuation functions" Ap q in the two-body 
Jastrow factor approximately removes its effect on one- 
body operators. Therefore, we assume initially that the 



coefficients v. 



(3) 

G 



0. We are unable to derive an analytic 



expression for the coefficients of the energy contri- 
bution e^ 2 - 1 and so, for the moment, we leave these aside. 
Constructing the full analytic approximation 



S G (x) 



(x), 



to an approximation to the coefficients of the functions 



«S(u) = 



TJ^qq'O? 2 + q' 2 ) + 2^W qk (k ■ k')(pk'-k>k'q' 
kk' 

Again, we are unable to derive expressions for the co- 
efficients corresponding to e' 2 '. We extract a two-body 
contribution from the constant contribution e^ 3 ' (see Ap- 
pendix ^), as u^, w — ^V(q)S(q' — q), for the electron 
interaction V. We replace the true interaction in this ex- 
pression with a pseudointeraction Vp S , which is generated 
for a given cut-off radius r c and reference eigenvalue e, as 
explained in Ref. V ps is used to generate the short- 
range two-body function u sr (Eq. 0), used in J sr . The 

(3) 

purpose of this modification of u qq , is to account for the 
presence of J sr in VP and is explained in Appendix ||]. 

From these contributions we construct the analytic 
predictor 

Sqq'(u) = V^(u) - ^p S (q)<5(q - q') • 

We notice that for periodic systems, this function is sep- 
arable in the points q of the first Brillouin zone (BZ). 
For each q in BZ, we may expand Sq+G.q+c as a func- 
tion of {uq+G,q+G'}, with no coupling to parameters 
w q '+H,q'+H' for q' q- This block diagonal form of the 
analytic predictor allows us to find the roots u 1 by solv- 
ing for each block (i.e., each q) individually. We do this 
using the Newton-Raphson method. 

A reliable initial guess for the Newton-Raphson 
method, rather than using u = 0, is the homogeneous 
If we regard the system as ho- 
: for G / 0, then the solution 



Solution Of iSqq' = 

mogeneous, i.e., (pg) 
is 



= %' - q) 



1 

AN 



8JW ps (q) 



(31) 



This 



'qq 



is used as the starting point for the Newton- 



Raphson iterations, and for all systems studied, this ini- 
tial guess produced convergent roots of 5 qq ' = 0. 



we see that the roots of Sg are trivially x = 0, as we 
would expect. 



2. Two-body Jastrow Factor 

According to Sec. 0, the two-body Jastrow correla- 
tion factor is divided into a short-range term J sr and an 
inhomogeneous term Jjh- We optimize the variational 
parameters u in Jjh. As for the one-body Jastrow, we 
expand the corresponding energy contribution e*- 1 ) which 
depends on Jjh alone (see Appendix ||]) . This leads us 



E. Numerical terms in the predictor 

To complete the construction of the predictor 
V^(a; a°) defined in Eq. |2^, for a given set of variational 
parameters a , we must determine some terms numeri- 
cally by fitting fluctuations in the energy of the system. 

We calculate V m (a.°) in Eq. by fitting the entire lo- 
cal energy, using Eq. [2^. Also, we use the fitting method 

to numerically determine the functions Vm\ot', cx°) given 
in Eq. |2^, by fitting the energy contribution e' 2 ' and each 
of its derivatives with respect to the parameters to be op- 
timized. Again, we note that is explicitly linear in 
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the parameters a and the approximation in Eq. ^8] is 
exact in this case. 

To do all this numerical work we use Monte Carlo sam- 
pling to determine the required expectation values: 

(AO m AO n ) 
(AEAO m ) 
(Ae^AO m ) 
(A(fg)A0 n ) 

for m, n ranging over the number of parameters N a in 
the set a. However, for simultaneous optimization of 
one- and two-body terms, we make some simplifications 
to reduce the computational workload. If we assume that 
the parameter a m can be varied independently of a n , this 
corresponds to assuming that (AO m AO„) = 0. 

In the expansion of the analytic local energy terms 
coming from Jih, outlined in Appendix |], we saw that 
these consisted of one- and two-body terms. However, 
the one-body terms contained Fourier coefficients of the 
average charge density (/3 q ), which are zero for q ^ G, a 
reciprocal latti ce vec tor. Also, the analytic form of the 
predictor (Sec. HID ) is separable in the k-points of the 
Therefore, we make the following 

AO n ): 



first Brillouin Zone 

approximations to the covariance matrix (AO 



(ACVG, q +G'A0 q ,+H', q '+H') = for q ^ q' ; 

(AOq+G,q+G' AOh) = for q ^ , 

for reciprocal lattice vectors G, G', H, H'. 

We regard one- and two-body optimization as indepen- 
dent for non-zero k-points in the first Brillouin Zone. If 
we also exclude the covariance terms between one- and 
two-body operators for q = 0, we find that this slows 
the convergence of the method for smaller systems. For 
larger systems, this exclusion prevents the convergence of 
the Newton- Raphson method at the first iteration of our 
method, thus halting the optimization process. However, 
we are not, in any way, confined to using the Newton- 
Raphson method to find the roots of the predictor, and 
other, more robust, root-finding methods might overcome 
this problem. 

Note that depends on the variational part (f> of the 
wave function which is being optimized. The variational 
components are Jib and Jih(q) for each q in BZ, where 
we define 



Jh(q) = exp 



GG' 



q+G,q+G'-Pq+G,q+G' 



so that J ih = n q J ih(q)- 

This greatly reduces the complexity of the predictor, 
without sacrificing convergence of the method for the sys- 
tems studied here. Ultimately, the predictor (a; a ) is 
itself only an approximation of the true function V m (ot), 
but by definition matches this function at the current 
values of the parameters being optimized. Therefore, 
approximations in the predictor affect only the rate of 
convergence of the method. 



The expressions for and its derivatives de^/dxc 
for the one-body Jastrow factor Jib are given in Ap- 
pendix [5]. Upon fitting these terms to the operators 
Og = Ap G , we may construct the function (a; a ). 
This defines the fitted terms Tq in the predictor (Sec. 
lHB| ). 

For a given q in BZ, we use the expressions for and 
the derivatives de^ / 9u q +G,q+G' > given in Appendix ^, 
corresponding to the two-body Jastrow factor (q) de- 

( 2) 

fined above. We construct the function v q4 _ G q+ c ( u q)> 
where u q = {it q+ G, q +G'; for G,G'}, and define the nU- 
merical predictor terms T q+ G, q +G' = V q+G q +c> wmcn 
contribute to the linear dependence of the predictor on 



IV. RESULTS 

We now apply this optimization method to diamond 
and rhombohedral graphite. We shall compare the cor- 
relation factors determined in both these systems, given 
the inhomogeneity and anisotropy of the electron charge 
density in graphite, relative to diamond. 

The convergence criterion of our iterative optimiza- 
tion: V m (a n ) = 0, requires the examination of possibly 
thousands of parameters, and it is difficult to visualize 
the overall convergence of the method. For illustrative 
purposes, we use the coefficients V m (a n ) to construct 
a single function V n . In the chosen Fourier basis, this 
amounts to reconstructing the real-space function V n 
from its Fourier coefficients. For one-body optimizations 



V n (r)=J2V G (x 



(32) 



G 



and for two-body optimizations 

U"(r,r')= (33) 
£ £ e- l ^ +G )- r U q+G , q+G ,(u")e l ^ G ')- r ' . 

q GG' 

Two-body functions, such as V n (r, r') and the Jastrow 
correlation function it(r, r'), are functions of six variables 
and so, extracting useful information from them is diffi- 
cult. For illustrative purposes, we indicate in Fig. [I] two 
points, A and B, in both the diamond and rhombohedral 
graphite structures, corresponding to high and low elec- 
tron charge density regions, respectively. A lies midway 
between two bonded carbon atoms and B lies midway 
between two layers of carbon atoms. We shall position 
the first electron at either A or B. The second electron 
shall be moved away from this position along one of the 
following line segments (indicated by heavy black lines 
in Fig. Q): A A' lying within a layer of carbon atoms; 
AA" perpendicular to the layers; BB 1 lying between two 
layers; and BB" perpendicular to the layers. 

By this means we may plot inhomogeneous two-body 
functions in terms of the relative separation of the elec- 
trons in the system. In particular, we may draw some 
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B 



B' 



(a) diamond 



(b) graphite 



FIG. 1: Crystal structure of (a) diamond and (b) rhom- 
bohedral graphite, illustrating stacked layers of hexagonally 
arranged carbon atoms for graphite and buckled layers for 
diamond. In both structures, the point A lies at the mid- 
point of a carbon-carbon bond, with the lines AA' and AA" 
extending within a layer and perpendicular to the layers, re- 
spectively. The point B lies midway between two layers (at 
a hexagonal interstitial point in diamond) and the lines BB' 
and BB" extend between the layers and perpendicular to the 
layers, respectively. 




12 3 4 
r (atomic units) 
(a) 



12 3 4 
r (atomic units) 
(b) 



FIG. 2: (a) The pseudo-interaction V^, s (solid line) and the 
Coulomb interaction V = 1/r (dashed line) versus electron 
separation r. (b) The short-range Jastrow function it S r, gen- 
erated from V ps , versus electron separation r, for angular mo- 
menta / = (solid line) and I — 1 (dashed line). V ps is 
generated using r c — 1.9 a.u. (indicated by vertical dotted 
line) and e = 0.2 hartree (see text). 



cell.1 



A. Removal of cusp 



conclusions about electron correlation in the system by 
examination of u. We may determine the isotropy of u 
by comparing plots of u with the first electron kept at 
the same point but the second electron moved in perpen- 
dicular directions, e.g., by comparing plots designated 
by A A' and AA" . The homogeneity of u may be seen 
by comparing plots of u with the second electron mov- 
ing in parallel directions from different positions of the 
first electron, e.g., by comparing AA' and BB' . Any dif- 
ferences between these plots of u are attributable to the 
inhomogeneity and anisotropy of the electron correlation 
factors in the systems studied. 

We use periodic boundaryj-.conditions (PBC) to ap- 
proximate the infinite crystalB The simulation cells con- 
sist of an N\ x N2 x N3 unit cell arrangement. The 
unit cells in each system are defined by the Bravais 
lattice basis vectors. For diamond, we use the ba- 
sis: {(a/2, a/2,0); (0, a/2, a/2); (a/2, 0, a/2)}, where 
a = 6.72 a.u., corresponding to a carbon bond-length 
of 2.91 a.u. For rhombohcdral graphite we use the basis: 
{(0,a,c); (-- \/3a/2,-a/2,c); (V^a/2, -a/2, c)}, where 
a = 2.68 a.u. is the bond-length within the layers, and 
c = 6.33 a.u. is the layer separation. 

We construct the Slater determinant D for both sys- 
tems using DFT-calculations in the local density approxi- 
mation (LDA).E3 The LDA orbitals were expanded using 
a linear combination of atomic orbitals comprising gaus- 
sians centred on each of the two carbon atoms in the unit 



In Sec. II A, we discussed the advantages of remov- 



ing the short-range cusp from the function u in order to 
improve its representability as a linear combination of 
smooth functions. Figure @(a) compares the pseudoint- 
eraction V^s, generated using a cut-off radius of r c = 1.9 
a.u. and energy eigenvalue e = 0.2 hartree (as discussed 
in Ref. |l|), with the Coulomb interaction V = e 2 /r. 
(In atomic units e 2 = 1). The pseudointeraction is used 
to generate a short-range Jastrow funtion u SI , which is 
shown in Fig. ||(b) for the relative angular momenta / = 1 
and I = corresponding to parallel spin and anti-parallel 
spin correlation, respectively. The short range Jastrow 
factor used in all subsequent calculations is that gener- 
ated with these particular values of r c and e. Subsequent 
figures in this paper, which involve u sr , represent anti- 
parallel spin correlation only. 

The cut-off required for a convergent Fourier expan- 
sion of a smooth cuspless function should be much less 
than that required for a function with a short-range cusp. 
Therefore, using a cuspless form greatly reduces the num- 
ber of terms required to represent the inhomogeneous 
form of the Jastrow factor accurately in Fourier space. 
We illustrate this point using a simple example. In Fig. 
^ we plot q 2 times the Fourier transform, for wave vector 
g, of the Yukawa-style homogeneous correlation factor, 



.4 



(1 



-r/F 



(34) 
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q 2 FT(u h -u sr ) 



4 6 
q (atomic units) 
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FIG. 3: Fourier transform times q 2 of (a) the homogeneous 
Jastrow factor «h (dashed line) and (b) the cuspless difference 
Mh — itsr (solid line) versus wave vector q. The homogeneous 
Jastrow factor parameter A — 1 a.u. (Eq. B4J) . 



which has been used by many authors to 
electron correlation in a variety of systems BBEIlL! 
We set A = 1 ajj_, (F is determined from A to satisfy the 
cusp conditionsEj and depends on the relative spin of the 
electrons). Also shown in Fig. |3| is q 2 times the Fourier 
transform of the cuspless difference Uh — u SI . We assume 
that, for large electron separations, electron correlation 
is approximately spin-independent. Therefore, for each 
g, we plot the mean of the parallel spin and anti-parallel 
spin values of the functions. 

In practice, we use the first zero of the Fourier trans- 
form of V ps as the Fourier space cut-off k c used in the 
definition of Jn, in Eq. [l4| For r c = 1.9 a.u. we use 
the cut-off k c — 2.185 (a.u.) -1 , beyond which the Fourier 
transform of the cuspless u function is approximately zero 
(Fig. |3|). Combining both the short-range and inhomo- 
geneous forms of the Jastrow factor using this scheme 
produces a form that is approximately independent of 
the cut-off, since decreasing r c increases the reciprocal 
space cut-off k c . 



B. Homogeneous Jastrow Factor 

For comparison with the inhomogeneous u functions 
determined in the following sections we use the Yukawa- 
style homogeneous function of Eq. [34| and construct 
a homogeneous trial wave function of the form ^ = 
J sl JhJibD. As with the inhomogeneous trial wave func- 
tion, we represent short-range correlation (i.e., the cusp) 
using J sr and represent the remaining correlation using a 
homogeneous Jastrow factor 



Jh = exp 



q<k a 




2 4 6 8 10 12 
r (atomic units) 

(a) diamond 



2 4 6 8 10 12 
r (atomic units) 



(b) graphite 



FIG. 4: Jastrow correlation factors versus electron separation 
r for (a) diamond and (b) rhombohedral graphite, (i) The 
function Uh (dotted line) as defined in Eq. B4J with optimized 
parameter A = 1.739 for diamond and 2.170 for graphite, 
(ii) The reconstructed function u of Eq. |35|, from a 3 x 3 x 3 
simulation region in both systems, for electron separations 
along the line segments AA' (solid line) and AA" (dashed 
line) from Fig. [j]. 



where we define 



u(q) 



[u h (r)-u SI (r)]e~^ r dr 



For the uniform electron gas, Bohm and PinesEl pre- 
dicted that the true function u should decay as l/w p r, 
at large separation r, where lo p is the plasma frequency. 
Rather than use this limiting value A — 1/lu p in Eq. 
|34| , it is common to treat A as a free parameter such 
that the energy is minimized. Using variational ,calcu- 
lations we can determine the optimal value of AE£l We 
optimize the one-body Jastrow factor Jib in Eq. [l5|, using 
our iterative method. The Jastrow factor J = J sr Jh Jib, 
combined with the Slater determinant D, is our best ap- 
proximation of the true many-body eigenstate using a 
homogeneous two-body Jastrow factor and is compara- 
ble with similaj waye functions optimized using variance 
minimization.LpEZIuj 

In Fig. |^ we plot Uh and compare it with the recon- 
structed function 



u(r) 



(r) + V 



q<k 



(35) 



for 3x3x3 unit cell simulations of diamond and rhom- 
bohedral graphite. Periodic boundary conditions and the 
anisotropy of the unit cell make this reconstructed form 
appear quite different to the original isotropic function 
Wh- 
in particular, for rhombohedral graphite the unit cell 
used is quite anisotropic, leading to marked differences in 
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(a) lxlxl 



(b) 1x1x2 




2 3 
iteration 



2 3 4 
iteration 



FIG. 5: Energy of diamond, in hartree/atom, for each iteration of the optimization process and various simulation cell sizes 
using the inhomogeneous Jastrow factor Jm (solid line): (a) 1 x 1 x 1 unit cell, with N = 8 electrons and N a = 140 variational 
wave function parameters; (b) 1 X 1 X 2, N = 16, N a = 196; (c) 2 x 2 x 2, N = 64, N a = 532; (d) 3 x 3 x 3, N = 216, N a = 1235. 
Also shown (dashed line) are energies during optimization of the one-body Jastrow factor Jib (with 84 variational parameters) 
in combination with the homogeneous Jastrow factor Jh. Averages for each iteration are calculated using 10 5 Monte Carlo 
samples. 



the reconstructed function along the perpendicular line 
segments AA' and AA" . (The homogeneity of it n is pre- 
served and so we only plot the function for point A, since 
all other points are equivalent.) Note that the Jastrow 
u function is defined up to an arbitrary constant, much 
like a potential, since this constant affects only the wave 
function normalization and contributes nothing to the 
description of correlation. Therefore, it is of no conse- 
quence that the functional form of uu in Eq. |34| appears 
shifted above the reconstructed forms compatible with 
PBC. This is due to the removal of the constant Fourier 
coefficient of the correlation function u(G = 0) from the 
expansion of the reconstructed-function in Eq. [35| We 
note that the cusp conditionals are maintained by all 
forms. 

We optimize the one-body Jastrow Jib using the 
method described in Sec. Ill . For the diamond simu- 



lations we used a Fourier space cut-off G c = 5.0 (a.u.)~ , 
giving 84 variational one-body parameters. This is more 
than enough for an accurate represe ntati on of one-body 
terms in the wave function (see Sec. VA). For graphite 
simulations we used G c = 3.1 (a.u.) , giving 30 varia- 
tional one-body parameters. 

The values of the electronic energy per atom for var- 
ious optimizations of this homogeneous trial wave func- 
tion are shown in Figs. [| and [| Each iteration involved 
averaging over 10 5 Monte Carlo samples. However, this 



amount of averaging was more than enough for an ac- 
curate implementation of the method. The optimization 
of the one-body terms in the 3x3x3 simulation of 
graphite (Fig. |J) involved only 2.5 x 10 4 samples per it- 
eration and is well converged. On average, the gain in 
energy following one-body optimization is approximately 
1.0 mhartree/atom for diamond and 2.5 mhartree/atom 
for graphite. We note that the necessity for a one-body 
correction is a consequence of the inhqmpgeneity of the 
electronic charge density in the systemJlJ'ES and it is not 
surprising that the gain in energy is larger for the more 
inhomogeneous system, graphite, than for diamond. 

The Slater determinant used in the diamond calcula- 
tions of Fig. ||is composed of single-particle orbitals gen- 
erated from an LDA calculation. The LDA orbitals are 
linear combinations of gaussian basis functions of s, p 
and d symmetry, using three decays of 0.24, 0.797, and 
2.65. The exchange-correlation functional used was of the 
Ceperley-AlderQ form. The cui^pff of the Fourier space 
expansion of the charge densitycJ in the LDA calculation 
was 64 Rydberg. 

The graphite calculations shown in Fig. ^ use a Slater 
determinant generated using a gaussian basis-set with 
s and p symmetry only. Four decays were used: 0.19, 
0.474, 1.183, and 2.95. The orbitals were generated frcwa, 
LDA calculations incorporating the Hedin-LundqvistE3 
exchange-correlation functional, and a cut-off of 36 Ry- 
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FIG. 6: Energy of rhombohedral graphite, in hartree/atom, for each iteration of the optimization process and various simulation 
cell sizes using the inhomogeneous Jastrow factor Jih (solid line): (a) 1 x 1 x 1 unit cell, with N = 8 electrons and N a ~ 162 
independent parameters for J ih ; (b) 1 x 1 x 2, N = 16, N a = 272; (c)2x2x2,iV = 64, N a = 1020; (d) 3 x 3 x 3, N = 216, 
N a = 3136. Also shown (dashed line) are energies during optimization of the one-body Jastrow factor Jib (with 30 variational 
parameters) in combination with the homogeneous Jastrow factor Jh. Averages for each iteration are calculated using 10 5 
Monte Carlo samples. (The optimization of Jib for the 3x3x3 simulation used only 2.5 x 10 4 samples per iteration.) 



dberg for the F ourie r space expansion of the charge den- 
sity. (See Sec. V_A for a discussion of basis-set conver- 
gence of the total energy in graphite.) 

Figure illustrates the convergence during optimiza- 
tion of the one-body function x( r ) in Eq. |l| (where x is the 
accumulation of all one-body terms from each of the Jas- 
trow factors) for the 3x3x3 diamond calculation. The 
optimized function is statistically well-determined, and 
the change from the initial one-body function, x 1 = X° 



of Eq. |Sj (as defined in Refs. [l?] and 18), is well-defined. 
This alteration of the one-body function may be com- 
pared to simi lar calculations performed using variance 
minimization.! 1 ^rl 18 ! 

The number of parameters for optimization could be 
greatly reduced through exploitation of the crystal point- 
group symmetry of the structures involved. However, it 
is worth noting that the optimization process preserves 
the natural symmetry of the system (within statistical 
error) without such measures, illustrating that for non- 
symmetric systems with large numbers of parameters, 
this optimization process should be quite robust. 

The one-body function V"(r), reconstructed from the 
coefficients associated with the local energy Vb at itera- 
tion n according to Eq. |3^, is shown in Fig. |^. Clearly, 
this function decreases in magnitude, indicating a de- 
crease in the magnitudes of the Euler-Lagrange deriva- 
tives. Beyond the first iteration, V n (r) is of the same 
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FIG. 7: The diamond one-body function x versus position 
r along the line segment AA' [Fig. [j](a)]. x" indicates the 
one-body function used at iteration n of the optimization of 
Jib in the presence of Jh in Fig. ||(d). 



order of magnitude as its associated standard error, and 
so is statistically insignificant. Therefore, the method has 
essentially converged after only one iteration. The noisi- 
est regions of V n (r) correspond to regions of low density 
where the Monte Carlo sampling is less frequent. 
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FIG. 8: The reconstructed one-body function V n , defined in 
Eq. [w], versus position r as for Fig. |?]. 
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FIG. 9: The Jastrow factor wrpa versus electron separation 
r along the line segments indicated in Fig. [j](b), for a 3 x 3 x 3 
simulation of rhombohedral graphite. 



C. Inhomogeneous RPA Jastrow Factor 

The analytic guess for the t wo- body predictor func- 
outlincd in Sec. 



tions S. 



qq'( U ), 



HID 2 



leads to an in- 
homogeneous generalization of the RPA equations. The 
solution to these equations is the function U1 = urpa, 
shown in Fig. ^ for the 3x3x3 graphite simulation. We 
notice some inhomogeneity in urpa at intermediate- and 
long-range electron separations. A more homogeneous 
and isotropic mrpa was found for diamond, as we would 
expect since diamond possesses a more uniform electron 
density than graphite. 

In Figs. H and [| the first point on all solid curves 
indicates the total energy per atom in each simulation 
calculated using urpa- In comparison with the energy 
calculated using the homogeneous Jastrow function u^, 
we see that the inhomogeneous RPA trial wave function 



is at best comparable in accuracy with the optimized 
trial wave function with homogeneous two-body Jastrow 
factor, and often less accurate_.These results are different 
from those of Gaudoin et alxJl for model systems: they 
find that their inhomogeneous generalization of the RPA 
produces wave functions that yield lower energies than 
the homogeneous form. 



D. Inhomogeneous Optimal Jastrow Factor 

We simultaneously optimized the parameters in both 
the one-body Jastrow factor Jib and the fully inhomoge- 
neous form of the two-body Jastrow factor Jjh, using our 
iterative method. The convergence of the total energy 
per atom for various simulations of diamond and rhom- 
bohedral graphite are shown in Figs. |5| and ^. (The Slater 
determinants used in com binat ion with the homogeneous 
Jastrow factor Jh in Sec. [VB are also used here.) Con- 
vergence of the total energy is achieved in approximately 
three iterations in all cases and is stable. This is remark- 
able given that the system sizes range from 8 electrons 
and 140 independent parameters to 216 electrons and 
3136 independent parameters. Also, we used the same 
amount of Monte Carlo sampling, viz., 10 5 samples per 
iteration, to determine the required expectation values 
for all simulations. In all cases the fully inhomogeneous 
form of u allows us to determine more accurate trial wave 
functions with substantially lower energies than the ho- 
mogeneous trial functions. In general, the gain in energy 
through using an inhomogeneous rather than a homoge- 
neous wave function is of the order of 5 mhartree/atom 
for both diamond and rhombohedral graphite. 

Figure [h] illustrates the rapid convergence of the two- 
body wave function parameters u to their optimal values 
during the largest graphite optimization (3x3x3) and 
is typical of all the optimizations performed in both dia- 
mond and graphite. Beyond the third iteration, no clear 
distinction exists between subsequent sets of parameters. 
The optimal Jastrow function u is significantly different 
from both the RPA function u 1 and the homogeneous 
form it/, of Eq. |34|. The proof that this optimization suc- 
ceeds in minimizing the energy expectation value may be 
seen in Fig. [ll] (which comes from the same graphite cal- 
culation as Fig. |o|). Here, we plot the iterative decay of 
the two-body function V n (r,r') reconstructed from the 
total energy coefficients determined at each iteration, ac- 
cording to Eq. |33|. This clearly indicates the reduction 
to zero (within statistical noise) of the derivatives in the 
Euler-Lagrange equations, thus solving the energy mini- 
mization problem. 

For the largest simulation cells studied (3x3x3 unit 
cell arrangement containing 216 electrons), we compare 
the optimal Jastrow correlation functions u of diamond 
and rhombohedral graphite. Figures |lj and [l3] show the 
function u plotted with respect to electron separation on 
various line segments in the corresponding crystal struc- 
tures, as already explained. 
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FIG. 



10: Graphite Jastrow correlation functions u" with 
respect to electron separation r on the line segment AA' (Fig. 
[l]) for each iteration n during optimization of for a 3 x 3 x 3 
simulation region [Fig. |^(d)] . Also shown is the reconstruction 
of the homogeneous function uh (heavy dotted line) defined 
in Eq. |35l with optimized parameter A — 2.170. 
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FIG. 11: The reconstructed two-body function V n , as de- 
fined in Eq. versus electron separation r on the line seg- 
ment AA' , for each iteration n during the optimization out- 
lined in Figs. [](d) and [io[ 



DISCUSSION 



A. Energies 



For a direct comparison of the calculated energies of 
diamond and rhombohcdral graphite, we should (a) make 
some corrections based on the trial wave functions used 
and (b) include finite size and zero-point phonon energy 
corrections for the expected energy of the real infinite 
solid. The corrected energies for diamond and graphite 
are listed in Table |. 

The Jastrow factors used in both solids are comparable 
in their variational freedom, the only significant differ- 
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FIG. 12: The optimized diamond Jastrow correlation factor 
u as a function of electron separation r along the line segments 
indicated, for a 3 x 3 x 3 simulation region. 
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FIG. 13: The optimized rhombohedral graphite Jastrow cor- 
relation factor u as a function of electron separation r along 
the line segments indicated, for a 3 x 3 x 3 simulation region. 



TABLE I: Energies and energy corrections of diamond and 
rhombohedral graphite (in hartree/atom). (a) Total energy, 
determined by VMC and our optimization method, for a 
3x3x3 unit cell simulation region, using wave functions of 
similar variational freedom (see text), (b) Finite size correc- 
tion equal to the energy difference between an LDA calcula- 
tion for the 3x3x3 simulation region and an LDA calculation 
using a fully-converged k-point set-Zsee text), (c) Correction 
for the zero-point phonon energy.E3 (d) Total energy includ- 
ing the corrections. The numbers in parentheses indicate the 
statistical error in the last digits of the corresponding energy. 



diamond 



graphite 



3x3x3 

finite size correction 
zero-point energy 
total 



-5.712 95(14) 
-0.008 99 
0.006 65 
-5.715 29(14) 



-5.712 91(14) 
-0.006 56 
0.006 10 
-5.713 41(14) 
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ence being the size of the cut-off used for xib- However, 
a VMC calculation for diamond, using the same cut-off 
as in graphite (G c = 3.1 (a.u.) -1 , corresponding to 25 
variational parameters), resulted in an increase in energy 
of only 1.0± 0.3 mhartree/atom for diamond. Therefore, 
the high Fourier coefficients of Xib contribute little to the 
total energy of the system. The diamond VMC energy 
for the 3x3x3 simulation quoted in Table | was deter- 
mined using G c — 3.1 (a.u.)^ 1 as the cut-off for one-body 
terms. 

In graphite, the exclusion of d symmetry from the basis 
set used to construct the LDA orbitals in the Slater de- 
terminant is energetically more im porta nt. In add ition to 
the calculations described in Sees. 1VB and [VP , we also 
performed calculations for graphite using a gaussian basis 
set with s, p and d symmetry, and three gaussian decays: 
0.22, 0.766, 2.67. For the 3 x 3 x 3 simulation, including 
d symmetry reduces the total VMC energy by 7.2 ± 0.3 
mhartree/atom and also reduces the variance in the to- 
tal energy by 16%. Use of the Ceperley- Alder exchange 
correlation functional to generate the single-particle or- 
bitals in graphite, rather than the Hedin-Lundqvist form, 
made no difference (within statistical error) to the VMC 
energies, and neither did an increase in the cut-off for the 
Fourier space expansion of the LDA charge density from 
36 to 64 Rydberg. The graphite VMC energy for the 
3x3x3 simulation listed in Table | was calculated using 
a trial wave function very similar to that used for the 
calculation of the corresponding diamond VMC energy. 
The cut-off for the one-body Jastrow factor was G c = 3.1 
(a.u.) -1 and the Slater determinant comprised LDA or- 
bitals obtained using: (i) a basis set with d symmetry (as 
outlined above); (ii) the Ceperley- Alder exchange corre- 
lation functional; and (iii) a 64 Rydberg cut-off for the 
LDA charge density expansion. 

We generate finite size corrections for the 3x3x3 unit 
cell simulations of diamond and graphite, by calculat- 
ing the difference in energy between an LDA calculation 
which uses k-points compatible with periodic boundary 
conditions of a 3 x 3 x 3 simulation region and-,an LDA cal- 
culation using a fully converged k-point setE Comparing 
the change in energy between a 2 x 2 x 2 calculation and a 
3x3x3 calculation in both diamond and graphite, using 
LDA and VMC, we see that the change in VMC energy 
is about 80% of the LDA energy change in diamond, and 
70% in rhombohedral graphite. Perhaps more accurate 
estimates of the energy of the infinite solid may be ob- 
tained by implementation of a model periodic Coulomb 
interaction developed recently. Tests of this approach 
have dramaticalhz. reduced finite-size effects in the inter- 
action energyO'ErEiJ 

For diamond, we estimated the finite size correction to 
be —8.99 mhartree/atom, using a converged LDA calcu- 
lation with 220 k-points in the irreducible Brillouin zone. 
For rhombohedral graphite, incorporating d symmetry 
in the basis set (as described above), and using 189 k- 
points in the LDA calculation, we found the finite size 
correction to be —6.56 mhartree/atom. We also include 



the calculated zero-point phonon energies of diamond 
and graphite- which are 6.65 and 6.10 mhartree/atom, 
respectively.Ea 

Adding all these corrections to the calculated VMC 
energies (Table |), we estimate the energies of the in- 
finite solids to be -5.71529 ± 0.00014 hartree/atom for 
diamond and —5.71341 ±0.00014 hartree/atom for rhom- 
bohedral graphite. This appears to indicate that rhom- 
bohedral graphite is less stable than diamond. However, 
given the approximation of using LDA finite size correc- 
tions, we might expect a systematic error of the order 
of 2 mhartree/atom in each of these results. This in- 
dicates that at the VMC level, the solids diamond and 
rhombohedral graphite have very similar total energies. 
We note that in the atomic pseudopotential used in the 
calculations presented here p and higher angular momen- 
tum scattering are all included in the local potential. It 
is possible that the use of a separate d pseudopotential 
might slightly affect the relative energies in both systems. 

In order to determine the cohesive energy of a solid, 
we should subtract the energy per atom of the solid from 
the energy of the isolated atom, E c = E a — E s . However, 
when using approximate eigenfunctions, a reasonable es- 
timate of the cohesive energy is obtainable only by sub- 
tracting the energies estimated using similar trial wave 
functions. VMC energies are available for the carbon 
atom where-the trial wave function is of the Jastrow- 
Slater form.El The orbitals of the Slater determinant are 
optimized using energy minimization and a sophisticated 
Jastrow factor is optimized using variance minimization, 
yielding a VMC energy for the atom of -5.4372 ± 0.0001 
mhartree. Using this energy, we find the cohesive ener- 
gies to be 0.2781 ± 0.0002 hartree/atom for diamond and 
0.2762 ± 0.0002 hartree/atom for graphite. We regard 
this atomic trial wave function to be close in form and 
accuracy to our solid trial wave function. However, to re- 
main consistent with the inclusion of d symmetry in the 
basis set of our LDA calculations, we could also refer to 
a multiconfiguration trial wave function for the carbon 
atom which includes d excitations. The VMC energy of 
the atom, using this wave function, is —5.45061 ±0.00002 
hartree^l leading to estimates of the cohesive energies of 
diamond and rhombohedral graphite of 0.2647 ± 0.0002 
hartree/atom and 0.2628 ± 0.0002 hartree/atom, respec- 
tively. The experimental values are 0.271 hartree/atom 
for diamond, and 0.272 hartree/atom for graphite£3 

Considering the significant gain in energy obtained us- 
ing an inhomogeneous Jastrow factor rather than a homo- 
geneous form, one might wonder how much the difference 
between VMC and DMC energies has been reduced. Us- 
ing Jastrow correlation factors which include one-body 
terms and homogeneous two-body termSj-sjmilar in for. 



to that used in Sec. IV B, for bulk carbonES and silico 



typically leads to VMC correlation energies that are ap- 
proximately 90% of the corresponding DMC correlation 
energies, i.e., 



E\MG — Ehf 
Edmc — Ehf 



0.90 
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where -EVmc is the energy calculated using VMC, -Bdmc 
is the energy calculated using DMC, and S H f is the 
Hartree-Fock energy. (-Erf has been estimated using; just 
the LDA Slater determinant as trial wave function.H) We 
have observed that including inhomogeneity in the cor- 
relation factor leads to a gain in energy of approximately 
10 mhartree/atom in the largest simulations of both dia- 
mond and graphite. If we estimate that the VMC correla- 
tion energy, determined using our optimized wave func- 
tion with homogeneous correlation terms, obtains 90% 
of -Edmc — £hf ) then using the optimized wave func- 
tion with inhomogeneous correlation terms obtains 96% 
of the DMC correlation energy. This would indicate that 
the energy difference between VMC and DMC for these 
solids has been decreased by 60%. 

It is important to emphasize that, whatever physical 
conclusions we would like to draw from these calcula- 
tions, our principal aim has been to optimize the trial 
wave function for a given Hamiltonian, such that the ex- 
pectation value of the total energy is minimized. This 
aim has been achieved for all the systems studied. 



B. Correlation Factors 

Diamond, with a relatively homogeneous and isotropic 
electron charge density, exhibits an approximately homo- 
geneous, isotropic Jastrow correlation function u (Fig. 
|i2| ). At large electron separations (beyond 6 a.u.) we 
see that the electron correlation factor u in diamond is 
well approximated by homogeneous and isotropic func- 
tions, as all the curves plotted are quite similar. Only 
slight deviations from homogeneity exist at short and 
intermediate electron separations. This inhomogeneity 
may be seen by comparing u plotted with its fixed coor- 
dinate at different points: A A' and A A" are quite similar 
at short range, but clearly distinct from BB' and BB" 
in the same region. This inhomogeneity may have sig- 
nificant effects on the short-range pair-correlation func- 
tions calculated for diamond-like systems using VMC 
methods. ESlSO 

Graphite, with clear regions of high electron charge 
density and well-defined regions of very low electron 
charge density between its layers, is a highly inhomo- 
geneous and anisotropic structure. This is borne out in 
Fig. [l3|, where the function u differs considerably in vari- 
ous regions and in various directions. At short-range, u is 
surprisingly homogeneous in comparison with diamond. 
At intermediate separations the function displays both 
inhomogeneous and anisotropic behaviour. Given that 
the layer separation in these simulations is 6.33 a.u., this 
indicates that inhomogeneous correlation between adja- 
cent layers in the system is not insignificant. This may 
prove important for van der Waals interactions between 
the layers in graphite. 

At long range, the anisotropy of the graphite correla- 
tion factor is clearly shown in Fig. 13, where the corre- 
lation factors for electron separation vectors r' — r lying 



parallel to the graphite planes (AA' and BB') are dis- 
tinctly different from those with the separation vector 
perpendicular to the planes (AA" and BB"). Inhomo- 
geneity (i.e. an explicit dependence on the position r of 
the first electron) is displayed at long range in the dif- 
ferences between the function along A A' and BB'. We 
might expect this, given that the line AA' lies within 
the graphite planes, where the charge density is concen- 
trated, whereas BB' lies in the very low charge density 
region between the planes [see Fig. [|(b)]. On the other 
hand, when the separation vector is oriented perpendic- 
ular to the graphite planes (A A" and BB"), inhomoge- 
neous effects are substantially smaller at long range. 



C. Computational Details 

In order to reduce the complexity of the physical 
results in Sec. |y|, some computational details of the 
method were not discussed. We present and discuss some 
of these details in this section. 

(1) This iterative method is trivially parallelizable. 
To determine the expectation values required to con- 
struct the predictor, Monte Carlo sampling may be per- 
formed independently on many workstations and the re- 
sults combined. To obtain total energies of the accu- 
racy presented in this paper requires ~ 10 5 Monte Carlo 
samples. However, for optimization of the wave function 
using our method, this amount of sampling is also suffi- 
cient for accurate estimations of the ex pectat ion values 
required to construct the predictor (Sec. HIE). 

The extra time required to accumulate the various con- 
tributions to the predictor represents less than 5% of the 
computational time needed to calculate each energy sam- 
ple. There are two reasons for this: Firstly, because we 
have expressed the variational components of the Jas- 
trow factor as linear combinations of the operators O n 
of Sec. 



Ill 



These O m need only be evaluated once, for 
each electron configuration, in order to construct both 
the Jastrow factor and the predictor. Secondly, and more 
importantly, since we must calculate the total energy as a 
sum of various contributions defined by the Hamiltonian 
of the system, all the energy contributions needed to con- 
struct the predictor are already available, either directly 
or by some simple manipulation. It is also important to 
note that, despite the complex form of the Jastrow fac- 
tor outlined in Sec. [□], the amount of computational time 
spent evaluating it is still less than or equal to that spent 
evaluating the Slater determinant, for a given electron 
configuration. 

The largest calculations presented in this paper were 
performed on a Beowulf cluster of fifteen 500 MHz dual- 
processor workstations. For a 5 iteration optimization 
using 10 5 samples per iteration, these calculations took 
about 50 hours on this cluster. However, perhaps half 
this amount of sampling would have produced compa- 
rable results [see (3) below]. The required memory for 
storing all the expectation values necessary for this cal- 
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FIG. 14: Total graphite energy in hartree/atom versus itera- 
tion number for a 2 x 2 x 2 simulation region. The number of 
Monte Carlo samples per iteration ranges fromlO 4 to5xl0 5 . 
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FIG. 15: The reconstructed two-body function V 5 , versus 
electron separation r along the line segment BB' in graphite, 
from the fifth iteration of each of the calculations in Fig. [l4|. 



this noise on the wave function accuracy is not clear. To 
reduce the computational cost of these calculations we 
would prefer to do the least amount of sampling necessary 
to produce the required results. 

Figure |TJ illustrates the effect of various amounts of 
sampling on the convergence of the total energy in the 
optimization of the 2x2x2 simulation of rhombohedral 
graphite. At 10 4 samples per iteration, the estimation 
of the required expectation values, during the first itera- 
tion of our method, is too crude to produce a convergent 
root of the predictor using the Newton- Raphson method. 
This leads to a wave function, used in the second itera- 
tion, with many noisy parameters which arc more diffi- 
cult to optimize, as the convergence of the energy shows. 
However, for sampling involving 2 x 10 4 samples per it- 
eration, or more, we see that the convergence of the total 
energy is identical (within statistical accuracy). Com- 
parison of the optimized wave function parameters also 
reveals only small differences, indicating that the only 
means of determining the true benefits of more sampling 
is by examination of the fitted coefficients V m (cy. n ) at 
each iteration n. 

In Fig. [l5| we see, from the reconstruction V 5 (r, r') of 
the fitted coefficients for two-body optimization at the 
fifth and final iteration of the method, that the most 
sampling (5 x 10 5 samples per iteration) reduces the mag- 
nitude of the Euler-Lagrange derivatives the most, indi- 
cating that these wave function parameters are the most 
accurate. However, for practical purposes, there is lit- 
tle distinction between the accuracy of the wave function 
once we increase the sampling beyond 2 x 10 4 samples 
per iteration. Thus, the optimizations presented in Sec. 
[fV| may have used more computational time than was 
strictly necessary. However, more testing is required to 
determine the minimum amount of sampling as a func- 
tion of system size. 



VI. CONCLUSION 



culation was approximately 20 MB. 

(2) The Newton- Raphson method (see Appendix |c]), 
while quadratically convergent near a root of a multidi- 
mensional system, does possess some convergence prob- 
lems far from the root. In the calculations presented 
in this paper, we found that below a certain minimum 
amount of sampling, the noise in the estimated expecta- 
tion values used to construct the predictor caused diver- 
gence of the solution to Eq. ||(] using the Newton- Raphson 
method. This problem might be solved through the use 
of a more robust root-finding method for the predictor 
function. 

(3) From the analysis in Appendix we see that noise 
in the iterative method comes from the finite sampling 
used to estimate those expectation values (listed in Sec. 
HIE) required to construct the predictor. The effect of 



We have developed a generalized form of electron cor- 
relation factor for trial many-body wave functions of elec- 
trons in periodic solids. This form allows us to represent 
fully inhomogeneous electron correlation in real physi- 
cal periodic systems. It is computationally efficient to 
evaluate, since the electron cusp, which we express as 
a homogeneous correlation factor, is separated from the 
fully inhomogeneous form. 

We have also developed a rapidly convergent iterative 
method for the optimization of all variational parame- 
ters in these wave functions, minimizing the total en- 
ergy of the given system. It uses the accurate techniques 
of quantum Monte Carlo sampling to achieve this op- 
timization and has allowed new insights into the form 
of many-electron correlation in systems with highly in- 
homogeneous charge densities. We estimate the differ- 
ence in energies calculated using the optimal inhomoge- 
neous two-body correlation factor and the optimal ho- 
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mogeneous correlation factor to be approximately 60% 
of the difference between the DMC and VMC energies 
obtained with homogeneous 2-body correlation terms. 

In diamond, the optimal correlation factor is approx- 
imately homogeneous and isotropic, with some inhomo- 
geneity at short- and intermediate-range electron sepa- 
rations. This is consistent with its comparatively homo- 
geneous and isotropic electron charge density. Graphite 
has an optimal correlation factor which is quite homoge- 
neous at short-range electron separation, but is signifi- 
cantly inhomogeneous and anisotropic at intermediate- 
and long-range electron separations, as one might ex- 
pect from its highly inhomogeneous and anisotropic elec- 
tron charge density. Nevertheless, it is remarkable that 
despite very large inhomogeneity in the electroH-,p. 
correlation function, found by previous authors Ji^t' 
the ideal inhomogeneous Jastrow two-body term, calcu- 
lated for the first time here for diamond and graphite, 
displays relatively small inhomogeneity. Whether this 
conclusion can be extended to other systems (e.g. in- 
volving strongly correlated d-electrons) remains an open 
question. 
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Therefore, we may express the derivative of (A) as 

l d{A) = (*\AO m \*) {*\A\9) (gjCyg) 
2da m (*|tf> 

= (AO m )-(A)(O m ) 

= (AAA0 m ), 

where AA(R) = A(R) - (A) and A0 ro (R) = O m (R) 
(O m ). 

APPENDIX B: LEAST-SQUARES FITTING 

Our task is to minimize the integral 



X 



by choosing the appropriate parameters Eq and {14}. 
We note first that, at the minimum, 

k 

where (E) is the expectation value of the total energy 
and (Ok) is the expectation value of the operator Ok- To 
fulfill the minimization, we must set all the remaining 
first derivatives of \ 2 to zero, i.e., 



dx 2 
dVi 



-2<*| \n-E -J2 v k o k 0, |tf ) = 



APPENDIX A: DERIVATIVE OF OBSERVABLES 
WITH RESPECT TO WAVE FUNCTION 
PARAMETERS 

The derivative of the expectation value of an observ- 
able A, with respect to a parameter a m G a of the pa- 
rameterized wave function ^(ot), may be written as 

d(A) (jfrjAjjg) 
da m ~ da m (*|*) 

We assume that A is independent of the parameters a. 
Also, for a time independent system, we may in general 
express $ as a real function. Differentiating, we find that 



8 



da r , 



mm = 2mo m m 



da r 



where we define the operator associated with variations 
of a m as O m = d/da m , with local value 



Om(R) 



i a*(R) 



d 



*(R) da m da r , 
for the many-body configuration R. 



ln*(R) , 



for each Upon substitution of the minimum value of 
Eq, this leads to 

(Ed) - (E)(Oi) =5^Vfc[<OfcOi) - (O k )(Ot)} . 



Since, for any operators a, b we may say that 

(ab) - (a)(b) = (AaAfc) 

where Aa = a — (a), then we are lead to the conclusion 
that the least squares fitting is equivalent to solving the 
linear system 



^V- fe (AO fe AO ; ) = (AEAOt 



for each I. 



APPENDIX C: NEWTON-RAPHSON METHOD 

A n inte gral part of the iterative procedure outlined in 
Sec. [II B is the determination of the parameters ot that 
solve the system in Eq. [30| The determination of the 
roots of any multidimensional function can be trouble- 
some. In all the optimizations presented in this paper, 
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the predictor function V^(a; a ) is a quadratic function 
of a, whose coefficients arc determined analytically, or 
by numerical fitting at the point a . We ignore the im- 
plicit dependence of on cc° for the purpose of finding 
a root, and solve the system V^a) — for all m. 

We use the Newton-Raphson methodB to determine 
the roots. This is an iterative method of improving suc- 
cessive guesses for a root of a function. It involves the 
computation of the function and its Jacobian ma- 
trix of derivatives with respect to a, at each guess. The 
iterations continue until convergence of the solution is 
achieved within a predefined tolerance. 

The success of the Newton-Raphson method for multi- 
dimensional systems relies heavily on the proximity of 
the initial guess to the root we seek. For this reason, 
the initial guess chosen is normally the parameter set of 
the p revious iteration of the procedure outlined in Sec. 
IIIB. To find the first set of parameters, by finding the 



roots of the analytic predictor S m , we require a good ini- 
tial guess of the root. Since the one-body parameters 
Xg are expected to be small by construction, an initial 
guess of zero for all G was found to be sufficient to pro- 
duce a convergent solution to the first application of the 
Newton-Raphson method. 

For the two-body problem, we rescale the variables 
w qq ' to improve the convergence-, of the root-finding 
method. According to the RPAjIa the long-range be- 
haviour of the u function should take the form u(r) = 
l/(jj p r, where the plasma frequency for a homogeneous 
system with electron charge density n is uj p = V '47m. 
The charge density n is determined in the simulation re- 
gion to be N/fl, for N the number of electrons in the 
simulation volume Q. Therefore, for small q in Fourier 
space, u behaves like 



e^ 2 ) defined in Sec. [II C. For Jib, we have that 
1 V? Jib 



-Y 

2^ Jib 
-^(V?lnJ lb + |VilnJ lb | 

i 



G 



+ \ Y XG-G'{G - G')-G'xg> Ap G + constant . 

GG' 

The constant terms are not required, so we ignore them. 

The energy contribution cannot be expanded an- 
alytically as a linear combination of the functions A/9q 
since ty' = iP/Jib is not explicitly a function of these 
coordinates. 



,(2) 



^ Jib ' *' 

G i 



However, is explicitly linear in the parameters %, 
with derivatives 

which are independent of X- 



APPENDIX E: TWO-BODY ENERGY 
CONTRIBUTIONS 



u(q) = 



4tt 



1 



flLU n 



'Nq 2 



This indicates large relative differences between values of 
w(q) for small q. Therefore, in inhomogeneous systems, 
it would be appropriate to rescale the parameters w qq ' by 
multiplying by A|q||q'|, thus rendering all the variables 
of the same order of magnitude as the plasma frequency. 
This leads to a less pathological numerical problem for 
the Newton-Raphson method. Appropriate scaling must 
also be applied to the predictor function , in Sec. III. 
For the initial analytic guess of the roots of S'qq/ we use 
the homogeneous solution for it qq ' outlined in Eq. |3l| . 



APPENDIX D: ONE-BODY ENERGY 
CONTRIBUTIONS 



The variational part of iff associated with one-body 
terms is the Jastrow factor Jib (Eq. |l5| ) with parameters 
X = {xg}- We expand the energy contributions and 



The variational part of \& associated with two-body 
energy contributions is the inhomogeneous Jastrow fac- 
tor Jjh of Eq. [l4|, with parameters u = {it qq '}- (For 
periodic systems we use only « q +G,q+G'-) The energy 



contributions of Sec. [II C are expanded here. The 
contribution is dependent only on the form of Jih and 
is expanded as 



1 \- Vf 

2^ J ih 



^(V 2 lnJ ih + |V, lnJ ih | 2 ) 



We find that 
-i]Tv 2 lnJ ih = 



qq' 



Y 



N -1 



N 



■Y u ™ 



<? 2 (Pq 



Ap* + constant 



We retain only the linear combination of the functions 
Jqq'- The constant terms we may ignore, and the one- 
body terms (linear combinations of Ap*) we assume are 
compensated by terms in the one-body Jastrow Jib- 
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If we ignore the removal of one-body terms from the 
function P qq / = (Ap q Ap*,)[^] and consider using the 
function A/9 q Ap q , instead, then we find that 



2 ^ 



V.lnJihl 2 



(El) 



"qk(k • k')w k - q /p k -_ k Ap q Ap q , 

qq' kk' 



where we have ignored constant and one-body terms 
in the expansion. The product p k '_ k A/9 q Ap q / contains 
both two- and three-body terms, since we may rewrite 
Pk'-k as A/5 k /_ k + (p k /_ k ). We intend here to remove 
two-body fluctuations and regard three-body fluctua- 
tions as mueh. less significant, so we retain only the two- 
body termcll (p k '_ k ) Ap q Ap*, from the charge fluctua- 
tion products in Eq. El, i.e., 



- 2 E E U ^ ■ k ')"k'q' (Pk'-k) A/9 q Ap q , + ■ • • . 
qq' kk' 

Now we make the assumption that removing one-body 
terms from this expression , i.e., replacing A/9 q A/9*, with 
P qq ' is a good approximation, and obtain the expression 



for v q q/(u) given in Sec. 



HID 2 



The contribution cannot be expanded analytically 
in the basis of fluctuation functions Pqq', However, it is 
clear that is linear in u, since, for \&' = ^//J;h, 



.(2) 



E u <; 

qq' 



E v ^ 



and has first derivatives 

9e (2) 



dUr 



which are independent of the parameters u. 

The final energy contribution we attempt to ex- 
press analytically in terms of the two-body fluctuation 
functions P qq '- If 'J' « D, the LDA Slater determi- 
nant, then the sum of contributions from the external po- 
tential and the kinetic energy term — (1/2) J^. V|\E , '/\E'' 
is a one-bodjz_.contribution defined by the Kohn-Sham 



Hamiltonianjfil since 



1V?D 



2 D 



+ wo] = E^ s - - v ^i)\ 



where ef s are the Kohn-Sham eigenvalues, Vh is the 
Hartree potential and V xc is the exchange and correla- 
tion potential. Therefore, the two-body contribution of 
these terms is approximately zero. We are left with the 
contribution of the electron-electron interaction V. 

For a two-body potential V(r, r'), we may expand the 
sum over electron pairs as 

X>(*i,ri) = jE^'E e ^ rieiq '^ 



for Fourier coefficients Vq q ' . In terms of the fluctuation 
functions P qq ' , this may be rewritten as 

J2v(r i ,r j ) = ^V^P^ 



E 



qq 

N - 1 
N 



constant 



(We assume that V(r, r') possesses exchange symmetry, 
so that Vqq' = V-q'-q.) Therefore, both one- and two- 
body fluctuations arise from a two-body potential in the 
"charge fluctuation" coordinate system. Note that the 
one-body fluctuations are expressible in terms of the 
Hartree potential, since 

q' 

If the two-body potential is homogeneous, i.e. 
V(r, r') = V(r — r'), then we may simplify the fluctua- 
tions since Vq q ' = V^S(q — q'), where the Fourier trans- 
form of the homogeneous function V(r) is 



V n 



V(r)e- 4qr dr 



for a system volume ft. If V(r) = V(— r) then V* = Vq. 
Therefore, 

E^i^-^i) = \ Y v * p ™ ( E2 ) 

q 

JV-1^„, . A , 
- — ^— Vq(pq) Ap q + constant . 



Again we ignore the one-body and constant terms in this 
context. 



APPENDIX F: CONSEQUENCES OF USING THE 
SHORT-RANGE JASTROW 

The short range Jastrow factor J sr is constructed from 
the pseudointeraction V" ps in the following way. For the 
isolated two-electron scattering problem, we may find an 
eigenstate ipo of the true two-electron Hamilatonian ho 
for a given energy eigenvalue e. Upon replacing the true 
Coulomb interaction V with a generated pseudointerac- 
tion V vs , we construct the modified Hamiltonian h ps with 
eigenstate ipps corresponding to e. We construct J sr such 
that 



^0 = Jsrlp 



ps 



qq' 



Now, in a many-electron environment, we know that 
using J sr allows for good approximation of short-range 
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correlations. We might imagine that for a many- 
electron system with Hamiltonian Ti. and many-electron 
trial wave function J sl ^ ps , the true energy eigenvalue 
may be well approximated by 



E 



where Ti ps = H- Y^i<jW{ r ij) ~ ^ps( r »i)]- Now > disguis- 
ing the true interaction V with V ps + ( V — V ps ) , and given 
the transferability of V ps over a wide range of energies, 
we see that 



^2V ps (nj) 



1 V? J S] 



E 



E 



'/si 



ps 



Dividing two-body correlation into short-range and in- 
homogeneous terms, we use a Jastrow factor of the form 
J sr Jih- The local energy determined using the trial wave 
function = J sr Jih^', where W = 1 J r /( J sr Jih), may be 
expanded as 



E 



1 Vf J s , 



"SI 



E 



Jsr */Jsr 



E y ( r «) 



f \ - Vj^a. \ - V»Jih v»^' 1 \ - vftg 
2 ^ J ih ^ J ih 2 ^ 1" 

+ ^Fext^) 

V^tw \ l^V 2 J ih v-ViJih Vi*' 

- aE^ + E^t^). 

For this reason, we use V ps in the expansion of the lo- 
cal energy for Jih to implicitly include the short-range 
Jastrow factor J sr . 



APPENDIX G: CONVERGENCE OF THE 
METHOD 

The convergence criterion V m (a) — 0, is numerically 
never exactly achieved. Given that the predictor con- 
tains some terms determined by statistical fitting, finite 
sampling errors exist, and this noise is passed on to the 
fitted parameters in from Eqs. E3 and E71 



We consider V n = {V m (a n )}, the fitted coefficients 
of the local energy at the n th iteration of the method. 
The method may be regarded as an iterative map M, 
such that the coefficients are determined via V„ + i = 
M(V„). There are two sources of noise in V„ + i: (i) 
noise inherited from V n , which produced the parameters 
<y. n+1 , which were used to construct the wave function 
1 &(a n+1 ), with which we evaluated the expectation values 
used to calculate V n +i; and (ii) finite sampling noise in 
the evaluation of the expectation values in Eqs. |2^ and 
|27| via Monte Carlo sampling. Therefore, we associate a 
set of variances er 2 = {cr 2 ^ for each m at step n}, arising 
from these two sources of noise, to the coefficients V„. 
The variance due to finite sampling alone, at each step 
n. is s 2 and we use the initial condition cr\ = s\ . This 
implies that the variance obeys the following iterative 
map: 



' n+l 



D n+1 



|A„ 



where A„ = VyM(V n ). For a convergent map M, we 
are guaranteed that |A„| 2 < 1 at convergence. |A„| is 
a measure of the convergence rate of the map M, with 
|A„| w implying fast convergence and |A n | ~ 1 imply- 
ing slow convergence. Note that if |A„| > 1 the map is 
divergent. 

Therefore, if we regard s„ rs s and A„ w A for all n, 
for constants s and A, then the converged value of the 
variance in the fitted coefficients is 



f 



1-lAl 2 



(Gl) 



We conclude that, given a convergent method M, the 
presence of statistical noise does not lead to successively 
more ill-determined parameters a, since their variance is 
also convergent. 

Note that Eq. Gl indicates that the variance in the 
fitted coefficients V m (a) is always greater than or equal 
to the variance estimated using finite sampling. However, 
our implementation of the method represented by the 
map M indicates that |A| <C 1, since we find that the 
majority of the coefficients V m (ct) ultimately end up with 
magnitudes approximately equal to their finite sampling 
errors, signifying that statistically they are zero. 

Therefore, the final conclusion to be drawn from Eq. 
Gl is that the accuracy of our optimization method de- 
pends ultimately on the finite sampling error. Therefore, 
increasing the computational workload, by increasing the 
amount of sampling, will result in more accurate opti- 
mizations of the wave function. This is demonstrated by 
the results in Sec. 
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